In [1]:
# import data
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPoint
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt
import folium
from folium.plugins import MarkerCluster
import geopandas as gpd
import numpy as np
import gurobipy as gp
from gurobipy import GRB
In [2]:
# read data
stations = gpd.read_file("Data/Fire_Stations_SBC.shp")
county = gpd.read_file("Data/tl_2019_06083_faces.shp")
In [3]:
# remove water part
county = county.loc[county["LWFLAG"] == 'L'].reset_index(drop=True)
county
Out[3]:
| TFID | STATEFP10 | COUNTYFP10 | TRACTCE10 | BLKGRPCE10 | BLOCKCE10 | SUFFIX1CE | ZCTA5CE10 | UACE10 | PUMACE10 | ... | METDIVFP | CNECTAFP | NECTAFP | NCTADVFP | LWFLAG | OFFSET | ATOTAL | INTPTLAT | INTPTLON | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 257492722 | 06 | 083 | 002211 | 2 | 2004 | None | 93454 | None | 08301 | ... | None | None | None | None | L | N | 220098 | +34.9497690 | -120.3695124 | POLYGON ((-120.37608 34.95466, -120.36430 34.9... |
| 1 | 219296658 | 06 | 083 | 003102 | 1 | 1258 | None | 93436 | None | 08302 | ... | None | None | None | None | L | N | 255 | +34.4938573 | -120.4754525 | POLYGON ((-120.47580 34.49391, -120.47569 34.4... |
| 2 | 219296666 | 06 | 083 | 003102 | 1 | 1255 | None | 93436 | None | 08302 | ... | None | None | None | None | L | N | 1493 | +34.4752880 | -120.4521852 | POLYGON ((-120.45261 34.47496, -120.45242 34.4... |
| 3 | 219296667 | 06 | 083 | 003102 | 1 | 1257 | None | 93436 | None | 08302 | ... | None | None | None | None | L | N | 113683 | +34.4706349 | -120.4566533 | POLYGON ((-120.45865 34.46769, -120.45850 34.4... |
| 4 | 219296675 | 06 | 083 | 003102 | 1 | 1263 | None | 93436 | None | 08302 | ... | None | None | None | None | L | N | 166285 | +34.4491085 | -120.4587685 | POLYGON ((-120.46231 34.45021, -120.46228 34.4... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 17274 | 219310979 | 06 | 083 | 001800 | 1 | 1089 | None | 93254 | None | 08302 | ... | None | None | None | None | L | N | 16549 | +34.9985013 | -119.8130570 | POLYGON ((-119.81407 34.99936, -119.81388 34.9... |
| 17275 | 219311368 | 06 | 083 | 001306 | 1 | 1014 | None | 93105 | 79282 | 08303 | ... | None | None | None | None | L | N | 149777 | +34.4230712 | -119.7456431 | POLYGON ((-119.74828 34.42247, -119.74817 34.4... |
| 17276 | 219313344 | 06 | 083 | 001906 | 4 | 4045 | None | 93105 | None | 08302 | ... | None | None | None | None | L | N | 175839 | +34.5812649 | -119.9848855 | POLYGON ((-119.98708 34.58219, -119.98700 34.5... |
| 17277 | 219313538 | 06 | 083 | 002906 | 2 | 2004 | None | 93117 | 79282 | 08303 | ... | None | None | None | None | L | N | 104502 | +34.4500980 | -119.8355449 | POLYGON ((-119.83790 34.44967, -119.83785 34.4... |
| 17278 | 219312610 | 06 | 083 | 002005 | 2 | 2001 | None | 93455 | 79417 | 08301 | ... | None | None | None | None | L | N | 115 | +34.8613367 | -120.4128687 | POLYGON ((-120.41296 34.86140, -120.41282 34.8... |
17279 rows × 45 columns
In [4]:
county_shape = county.unary_union
In [5]:
# plot data
A = stations.plot(color="red")
gpd.GeoSeries(county_shape).plot(ax=A)
Out[5]:
<Axes: >
In [6]:
stations.crs
Out[6]:
<Projected CRS: EPSG:2229> Name: NAD83 / California zone 5 (ftUS) Axis Info [cartesian]: - X[east]: Easting (US survey foot) - Y[north]: Northing (US survey foot) Area of Use: - name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura. - bounds: (-121.42, 32.76, -114.12, 35.81) Coordinate Operation: - name: SPCS83 California zone 5 (US Survey feet) - method: Lambert Conic Conformal (2SP) Datum: North American Datum 1983 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
In [7]:
county.crs
Out[7]:
<Geographic 2D CRS: EPSG:4269> Name: NAD83 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands. - bounds: (167.65, 14.92, -40.73, 86.45) Datum: North American Datum 1983 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
In [8]:
# convert crs
county_shape = gpd.GeoDataFrame(geometry=[county.to_crs(stations.crs).unary_union], crs=stations.crs)
county_shape
Out[8]:
| geometry | |
|---|---|
| 0 | MULTIPOLYGON (((5776084.957 2033858.920, 57760... |
In [9]:
# Now, let's plot them together.
A = stations.plot(color="red", zorder=100, markersize=10)
county_shape.plot(ax=A, facecolor='None', edgecolor="grey")
# Review for SCALEBAR!
scalebar = ScaleBar(dx=1, scale_formatter=lambda data, unit: f'{data},000 ft', length_fraction=0.25)
A.add_artist(scalebar)
Out[9]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x21eac736410>
In [10]:
import folium
from folium.plugins import MousePosition
import geopandas as gpd
county_shape_4326 = county_shape.to_crs(4326)
# Calculate the centroid of the polygon
centroid = county_shape_4326.geometry.centroid.values[0].coords[0]
# Create a Folium map centered at the centroid of the polygon
m = folium.Map(location=[centroid[1], centroid[0]], zoom_start=9)
# Add the polygon to the map
folium.GeoJson(county_shape_4326, style_function=lambda x: { 'fillColor': 'transparent'}).add_to(m)
folium.GeoJson(stations.to_crs(4326)).add_to(m)
# Add MousePosition control to display cursor coordinates
#### Use the cursor coordinates to define your Area of Interest! ####
MousePosition().add_to(m)
# Display the map
m
C:\Users\18428\AppData\Local\Temp\ipykernel_13564\4168181903.py:9: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. centroid = county_shape_4326.geometry.centroid.values[0].coords[0]
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
# changing boundary
max_x, min_x = -120.39153, -120.49805
max_y, min_y = 34.98936, 34.84762
# Create the rectangle geometry, which forms are study region boundary.
boundary = Polygon([(min_x, min_y), (min_x, max_y), (max_x, max_y), (max_x, min_y), (min_x, min_y)])
In [12]:
folium.GeoJson(boundary.__geo_interface__, style_function=lambda x: {'color': 'red', 'fillColor': 'transparent'}).add_to(m)
m
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [13]:
boundary = gpd.GeoSeries(boundary, crs=4326).to_crs(stations.crs)[0]
boundary.area / 2.788e+7
Out[13]:
59.08324173317982
In [14]:
stations.within(boundary)
Out[14]:
0 False
1 False
2 False
3 False
4 False
...
203 False
204 False
205 False
206 False
207 False
Length: 208, dtype: bool
In [15]:
study_stations = stations.loc[stations.within(boundary)].reset_index(drop=True)
len(study_stations)
Out[15]:
9
In [16]:
boundary_county = county_shape.intersection(boundary)# the output will be assigned to boundary_county variable
boundary_county
Out[16]:
0 MULTIPOLYGON (((5817885.812 2191416.585, 58183... dtype: geometry
In [36]:
COVERAGE_DISTANCE = 2 * 5280 # because 1*5280 feet is 1 mile
# Draw buffer from each fire station in study_stations GeoDataFrame
# Buffer distance is COVERAGE_DISTANCE we've defined
stations_coverage = study_stations.buffer(COVERAGE_DISTANCE)
# Let's plot the county polgyon within our study region, and
# stations_coverage buffers
A = stations_coverage.plot(facecolor="None")
boundary_county.plot(ax=A, zorder=-1)
plt.axis("off")
# How much area is covered?
# To get to know, firstly take unary_union of the total station_coverage
total_stations_coverage = stations_coverage.unary_union
# Then take the intersection with boundary_county, to exclude water area
total_stations_coverage = total_stations_coverage.intersection(boundary_county)
# Then divide the area of the updated total_station_coverage with the boundary_county.area
coverage = total_stations_coverage.area / boundary_county.area
# Multiply 100 to get the % value
coverage = coverage * 100
# Check out how much percentage of study region is covered by current fire stations.
# Do they agree with each other? the coverage % and plot?
print(f'{round(coverage.values[0],2)} % covered')
76.7 % covered
In [37]:
redundant_coverage = 0
# This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
for i in range(len(stations_coverage)-1):
# Take one station coverage buffer
cover_i = stations_coverage[i]
for j in range(i+1, len(stations_coverage)):
# Take another station coverage buffer which hasn't paired with i
cover_j = stations_coverage[j]
print("\t combination: ", i, j)
# With the two stations combination, calculate the intersection area.
# and sum up!
redundant_coverage += cover_i.intersection(cover_j).area
print() # To add blank line
print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
combination: 0 1 combination: 0 2 combination: 0 3 combination: 0 4 combination: 0 5 combination: 0 6 combination: 0 7 combination: 0 8 combination: 1 2 combination: 1 3 combination: 1 4 combination: 1 5 combination: 1 6 combination: 1 7 combination: 1 8 combination: 2 3 combination: 2 4 combination: 2 5 combination: 2 6 combination: 2 7 combination: 2 8 combination: 3 4 combination: 3 5 combination: 3 6 combination: 3 7 combination: 3 8 combination: 4 5 combination: 4 6 combination: 4 7 combination: 4 8 combination: 5 6 combination: 5 7 combination: 5 8 combination: 6 7 combination: 6 8 combination: 7 8 Redundant Coverage: 75.52173527741837 square miles
In [38]:
def calculate_redundant_coverage(coverage_buffers):
redundant_coverage = 0
# This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
for i in range(len(coverage_buffers)-1):
# Take one station coverage buffer
cover_i = coverage_buffers[i]
for j in range(i+1, len(coverage_buffers)):
# Take another station coverage buffer which hasn't paired with i
cover_j = coverage_buffers[j]
# With the two stations combination, calculate the intersection area.
# and sum up!
redundant_coverage += cover_i.intersection(cover_j).area
print() # To add blank line
print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
return redundant_coverage / 2.788e+7
In [39]:
current_rc = calculate_redundant_coverage(stations_coverage)
Redundant Coverage: 75.52173527741837 square miles
In [40]:
boundary.bounds
Out[40]:
(5812247.644516614, 2139386.8749094545, 5845424.880169843, 2191733.8904521973)
In [41]:
minx, miny, maxx, maxy = boundary.bounds
print(minx, miny, maxx, maxy)
5812247.644516614 2139386.8749094545 5845424.880169843 2191733.8904521973
In [42]:
interval = 5280/3 # default: 5280/5 feet (1 mile / 5 = 0.2 mile)
# Create arrays of x and y coordinates using np.arange
x_coords = np.arange(minx, maxx, interval)
y_coords = np.arange(miny, maxy, interval)
# Create a list to store the points
fishnet_points = []
# Generate points for the fishnet
for y in y_coords:
for x in x_coords:
fishnet_points.append((x, y))
# Print the number of points generated
print("Number of points in the fishnet:", len(fishnet_points))
fishnet_points = gpd.GeoSeries([Point(pt_cd) for pt_cd in fishnet_points])
fishnet_points.plot(figsize=(15,5))
Number of points in the fishnet: 570
Out[42]:
<Axes: >
In [43]:
fishnet_points = fishnet_points.loc[fishnet_points.intersects(boundary_county[0])]
fishnet_points = fishnet_points.loc[fishnet_points.intersects(total_stations_coverage[0])].reset_index(drop=True)
fishnet_points.plot(figsize=(15,5))
Out[43]:
<Axes: >
In [44]:
A = fishnet_points.plot(figsize=(15,5), color="black", markersize=10)
study_stations.plot(ax=A, color="blue", marker="*", markersize=400)
Out[44]:
<Axes: >
In [45]:
# optimize the problem
num_demands = len(fishnet_points)
num_facilities = len(study_stations)
D_ij = [[np.sqrt((fishnet_points.x[i] - study_stations.geometry.x[j])**2 + (fishnet_points.y[i] - study_stations.geometry.y[j])**2)
for j in range(num_facilities)] for i in range(num_demands)]
N_i = dict(zip(range(num_demands), [[] for _ in range(num_demands)]))
for i in range(num_demands):
for j in range(num_facilities):
if D_ij[i][j] <= COVERAGE_DISTANCE:
N_i[i].append(j)
# Create a new model
model = gp.Model("LSCP")
# Decision variables
# x[i] = 1 if facility i is selected, 0 otherwise
x = model.addVars(num_facilities, vtype=GRB.BINARY, name="x")
# Objective function: minimize the number of selected facilities
model.setObjective(gp.quicksum(x[j] for j in range(num_facilities)), GRB.MINIMIZE)
# Constraints
# Each demand point must be covered by at least one selected facility
for i in range(num_demands):
model.addConstr(gp.quicksum(x[j] for j in N_i[i]) >= 1, name=f"cover_demand_point_{j}")
# Optimize the model
model.optimize()
Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (win64 - Windows 11.0 (22631.2)) CPU model: 11th Gen Intel(R) Core(TM) i7-11370H @ 3.30GHz, instruction set [SSE2|AVX|AVX2|AVX512] Thread count: 4 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 383 rows, 9 columns and 794 nonzeros Model fingerprint: 0xbcb456a9 Variable types: 0 continuous, 9 integer (9 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+00, 1e+00] Bounds range [1e+00, 1e+00] RHS range [1e+00, 1e+00] Found heuristic solution: objective 8.0000000 Presolve removed 383 rows and 9 columns Presolve time: 0.00s Presolve: All rows and columns removed Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units) Thread count was 1 (of 8 available processors) Solution count 1: 8 Optimal solution found (tolerance 1.00e-04) Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0000%
In [46]:
selected_bool = [x[i].x > 0.5 for i in range(num_facilities)]
selected_facilities = study_stations.loc[selected_bool].reset_index(drop=True)
unselected_facilities = study_stations.loc[~np.array(selected_bool)].reset_index(drop=True)
print("Objective value: ", model.objVal, sum(selected_bool))
print("Current running stations: ", len(study_stations))
Objective value: 8.0 8 Current running stations: 9
In [47]:
A = gpd.GeoSeries(MultiPoint(fishnet_points)).plot(figsize=(15,5), color="grey", markersize=5)
unselected_facilities.plot(ax=A, color="blue", legend=True, label="Unselected", marker="*", markersize=300, zorder=10)
selected_facilities.plot(ax=A, color="red", marker="*", markersize=300, legend=True, label="Selected")
selected_facilities.buffer(COVERAGE_DISTANCE).plot(ax=A, edgecolor="red", facecolor="None", linestyle="dashed")
plt.legend()
plt.axis("off")
Out[47]:
(5815990.983205412, 5853511.824250226, 2132101.015651713, 2202523.549244632)
In [48]:
# Convert GeoPandas GeoSeries to GeoDataFrame
fishnet_gdf = gpd.GeoDataFrame(geometry=fishnet_points, crs=stations.crs).to_crs(4326)
# Create a base map
m = folium.Map(location=[fishnet_gdf.unary_union.centroid.y, fishnet_gdf.unary_union.centroid.x], zoom_start=11)
# Plot fishnet
for idx, row in fishnet_gdf.iterrows():
folium.CircleMarker(location=[row.geometry.y, row.geometry.x],
radius=1,
color='grey', fillcolor="grey").add_to(m)
# Plot study stations
for idx, row in unselected_facilities.to_crs(4326).iterrows():
folium.Marker(location=[row.geometry.y, row.geometry.x],
icon=folium.Icon(color='blue'),
popup=row['Label']).add_to(m)
# Plot selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
folium.Marker(location=[row.geometry.y, row.geometry.x],
icon=folium.Icon(color='red', icon='star'),
popup=row['Label']).add_to(m)
# Plot buffer around selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
folium.vector_layers.Circle(location=[row.geometry.y, row.geometry.x],
radius=COVERAGE_DISTANCE/3.281, # feet to meter conversion
color='red',
fill=False,
dash_array='10').add_to(m)
# Display the map
m
Out[48]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [49]:
optimized_rc = calculate_redundant_coverage(selected_facilities.buffer(COVERAGE_DISTANCE))
Redundant Coverage: 51.62504510669677 square miles
In [50]:
spatial_efficiency = model.objVal / len(study_stations) * 100
print("Spatial Efficienty is ", round(spatial_efficiency, 2), " %")
Spatial Efficienty is 88.89 %
In [51]:
print("Coverage Distance (mi): " , COVERAGE_DISTANCE) # Double Check this value before putting into table **
print("Current Stations #: " ,len(study_stations))
print("Optimized Stations #: ", len(selected_facilities))
print("Current Redundant Coverage (sqmi): ", current_rc)
print("Optimized Redundant Coverage (sqmi): ", optimized_rc)
print("Spatial Efficiency (%): ", spatial_efficiency)
Coverage Distance (mi): 10560 Current Stations #: 9 Optimized Stations #: 8 Current Redundant Coverage (sqmi): 75.52173527741837 Optimized Redundant Coverage (sqmi): 51.62504510669677 Spatial Efficiency (%): 88.88888888888889
In [ ]: